home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / swagg_m.zip / MATH.SWG / 0085_Gaussian Distribution.pas < prev    next >
Pascal/Delphi Source File  |  1995-03-03  |  2KB  |  142 lines
  1. {
  2. >Now I need a fast way of generating random numbers in Gaussian Distribution
  3.  
  4. This is my code for gaussian and uniform distribution.
  5. After the unit there is a graphic test program.
  6.  
  7. From: randyd@alpha2.csd.uwm.edu (Randall Elton Ding)
  8. }
  9.  
  10. unit rndgauss;
  11.  
  12.  
  13. interface
  14.  
  15.  
  16. function rnd: double;                  { returns uniform [0..1] }
  17.  
  18. function gauss(a,d: double): double;   { a is mean, d is std deviation }
  19.  
  20.  
  21. implementation
  22.  
  23.  
  24. function rnd: double;
  25.   const
  26.     bias = 1023;
  27.  
  28.   var
  29.     data: record
  30.             b: byte;
  31.             d: double;
  32.           end;
  33.     x: array[0..8] of byte absolute data;
  34.     e,i,j: word;
  35.  
  36.   begin
  37.     for i:= 0 to 7 do x[i]:= lo(random(256));
  38.     e:= bias;
  39.     repeat
  40.       j:= 0;
  41.       for i:= 0 to 7 do begin
  42.         j:= (x[i] shl 1) + hi(j);
  43.         x[i]:= lo(j);
  44.       end;
  45.       e:= e-1;
  46.       if (bias-e) mod 8 = 0 then x[0]:= lo(random(256));
  47.     until (x[7] and $10) = $10;
  48.     x[7]:= (x[7] and $0F) or lo(e shl 4);
  49.     x[8]:= lo(e shr 4);
  50.     rnd:= data.d;
  51.   end;
  52.  
  53.  
  54.  
  55. function gauss(a,d: double): double;
  56.   const
  57.     t: double = 0;
  58.  
  59.   var
  60.     v1,v2,r: double;
  61.  
  62.   begin
  63.     if t=0 then begin
  64.       repeat
  65.         v1:= 2*rnd-1;
  66.         v2:= 2*rnd-1;
  67.         r:= v1*v1+v2*v2
  68.       until r<1;
  69.       r:= sqrt((-2*ln(r))/r);
  70.       t:= v2*r;
  71.       gauss:= a+v1*r*d;
  72.     end
  73.     else begin
  74.       gauss:= a+t*d;
  75.       t:= 0;
  76.     end;
  77.   end;
  78.  
  79.  
  80.  
  81. begin
  82. end.
  83.  
  84.  
  85. {---------------------- cut ---------------------------------}
  86.  
  87.  
  88. program gaussiantest;
  89.  
  90. uses crt,graph,rndgauss;
  91.  
  92. const
  93.   bgipath = 'c:\bp\bgi';
  94.   largestx = 999;
  95.  
  96.  
  97. procedure testplot;
  98.   var
  99.     htarry: array[0..largestx] of integer;
  100.     x,y,w,h,m,v: word;
  101.  
  102.   begin
  103.     fillchar(htarry,sizeof(htarry),#0);
  104.     w:= getmaxx+1;
  105.     h:= getmaxy;
  106.     m:= getmaxx div 2;
  107.     v:= getmaxx div 8;
  108.     while not keypressed do begin
  109.       x:= trunc(gauss(m,v));
  110.       if x<=largestx then begin
  111.         y:= htarry[x];
  112.         if y<=h then begin
  113.           putpixel(x,h-y,white);
  114.           htarry[x]:= y+1;
  115.         end;
  116.       end;
  117.     end;
  118.   end;
  119.  
  120.  
  121.  
  122. procedure initbgi;
  123.   var errcode,grmode,grdriver: integer;
  124.   begin
  125.     grdriver:= detect;
  126.     initgraph (grdriver,grmode,bgipath);
  127.     errcode:= graphresult;
  128.     if errcode <> grok then begin
  129.       writeln ('Graphics error: ',grapherrormsg (errcode));
  130.       halt (1);
  131.     end;
  132.   end;
  133.  
  134.  
  135.  
  136. begin
  137.   initbgi;
  138.   testplot;
  139.   closegraph;
  140. end.
  141.  
  142.